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We study the stability properties of the Kidder-Scheel-Teukolsky (KST) many-parameter formu- 
lation of Einstein's equations for weak gravitational waves on flat space-time from a continuum 
and numerical point of view. At the continuum, performing a linearized analysis of the equations 
around flat spacetime, it turns out that they have, essentially, no non-principal terms. As a con- 
sequence, in the weak field limit the stability properties of this formulation depend only on the 
level of hyperbolicity of the system. At the discrete level we present some simple one-dimensional 
simulations using the KST family. The goal is to analyze the type of instabilities that appear as 
one changes parameter values in the formulation. Lessons learnt in this analysis can be applied in 
"^j", other formulations with similar properties. 
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>Y I. INTRODUCTION 

Numerical simulations of the Einstein equations for situations of interest in the binary black hole problem do not 
run forever. The codes either stop due to the development of floating point overflows, or even if they do not crash, 
— i ' they produce answers that after a while are clearly incorrect. It is usually very difficult to pinpoint a clear reason 
why a code stops working. Recently, Kidder, Scheel and Teukolsky (KST) |Q introduced a twelve-parameter family 
i of evolution equations for general relativity. Performing an empirical parameter study within a certain two-parameter 
subfamily, they were able to evolve single black hole spacetimes for over 1000 M, where M is the mass of the black hole, 
something that had been very difficult to achieve in the past. It is of interest to try to understand what makes some 
of the parameter choices better than others, in particular given that a twelve dimensional parameter study appears 
, prohibitive at present. The intention of this paper is to take some steps in this direction. We will first perform a 
linearized analysis of the KST equations in the continuum, by considering small perturbations around flat space-time. 
C\| , We will observe that the stability of flat space-time is entirely characterized by the level of hyperbolicity of the system. 

Since the latter is controlled by the parameters of the family, this provides a first analytic guidance as to which values 
to choose. Unfortunately, the result is somewhat weak, since it just points to an obvious fact: formulations with 
QT ! higher level of hyperbolicity work better. 

In the second part of the paper we perform a set of simple numerical tests. We consider spacetimes where all 
(si). variables depend on one spatial coordinate, which we consider compactified for simplicity, and time. We are able 
to exhibit explicitly the various types of instabilities that arise in the system. Some of the results are surprising. 
For the situation where the system is weakly hyperbolic, the code is strictly non-convergent, but it might appear to 
converge for a significant range of resolutions. We will see that the addition of dissipation does not fix these problems, 
I but actually can exacerbate them. It is often the case in numerical relativity that discretization schemes that are 
convergent for strongly hyperbolic equations are applied to weakly hyperbolic formulations. The examples of this 
section will teach us how dangerous such a practice is and confirm the analytic results of reference . This part of 
the paper is also instructive in that the KST system has only been evolved with pseudo-spectral methods. We use 
ordinary integration via the method of lines. 

The plan of this paper is as follows. In the next section we will discuss several notions of stability that are present 
in the literature, mostly to clarify the notation. In section III we discuss the stability of the KST equations in the 
continuum under linearized perturbations. In section IV we discuss the numerical simulations. 
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II. DIFFERENT DEFINITIONS OF STABILITY 

The term stability is used in numerical relativity in several different ways. We therefore wanted to make the notation 
clear at least in what refers to this paper. Sometimes the notion of stability is used in a purely analytic context, 
while some other times it is used in a purely numerical one. Within analytical contexts, there are cases in which it is 
used to mean well posedness, as in the book of Kreiss and Lorenz jjj. In such a context well posedness means that 
the norm of the solution at a fixed time can be bounded by the norm of the initial data, with a bound that is valid 
for all initial data. In other cases it is intended to measure the growth of perturbations of a certain solution within a 
formulation of Einstein's equations, without special interest in whether the equations are well posed or not. 

At the numerical level, a scheme is sometimes defined as stable if it satisfies a discrete version of well posedness. 
This is the sense in which stability (plus consistency) is equivalent to convergence via the Lax theorem Q] . Examples 
of this kind of instability are present in the Euler scheme, schemes with Courant factor that are too large, or other 
situations where the amplification factor (or its generalization) is bigger than one. Finally, the term stability has also 
historically been loosely used to mean that a simulation runs for a certain amount of time before crashing, or that 
the errors remain reasonably bounded as a function of time for a certain amount of time. 

At the continuum, part of the problem we want to look at can then be stated as follows. Suppose there is a certain 
solution uq of Einstein's equations which is bounded for all times. One can prescribe the initial data /o that uniquely 
determines, through the use of the evolution equations, that solution. However, if one gives initial data slightly 
different from /o, the corresponding solution might grow without bound or even blow up in a finite amount of time. 
This can happen either because this is the physical situation, or because the new solution u\ does not satisfy the 
constraints, or does satisfy them but it is in a gauge that is becoming ill defined. Numerical simulations may blow up 
even when one tries to model physically stable spacetimes. This is the case for evolutions of Schwarzschild or Kerr 
black holes (which arc known to be at least linearly stable with respect to physical perturbations.) 

The time growth of solutions of a given set of partial differential equations can be classified according to whether 
or not it can be bounded by a time function that does not depend on the initial data. This leads to the concept of 
well posedness. Namely, a system of partial differential equations is said to be well posed at a solution Uo(t) if there 
is a T and some norm such that 

\\u(t)\\ < F(T)||u(0)|| forO<t<T (1) 

for all solutions u(t) with u(0) sufficiently close to Uo(0)- For linear equations the bound can be tightened to an 
exponential, but for nonlinear systems F(T) might be any other function that can even diverge at finite T. The 
strength of equation (Q) characterizes the property that the bound is independent of the details of the initial data; 
in particular, of their frequency components. This is important in numerical simulations, where more (and higher) 
frequency components appear when resolution is increased. Well posedness is a necessary condition in order to have 
long term, convergent evolutions, but it is not sufficient. Indeed, the function F(T) can grow quickly in time. In 
order to control F(T) one has to go beyond well posedness and study non-principal terms. We will now apply these 
ideas to the KST equations. 

III. THE KST FAMILY OF EQUATIONS AND THEIR LINEAR STABILITY 

A. The formulation 

Starting from the Arnowitt, Deser and Misner (ADM) equations, KST derive a family of strongly hyperbolic 
first-order evolution equations for the three- metric (<?y ), the extrinsic curvature (JQj ), and the spatial derivatives of 
the three-metric (dkij = dk9ij), which generalizes previous well posed formulations ||, [(J. A priori prescribing the 
densitized lapse exp (Q) as a function of spacetime, the lapse N is given by 

N = g°eQ, (2) 

where g is the determinant of the three- metric and a a parameter. The shift vector, N l , is also assumed to be 
prescribed as a function of spacetime. We define the derivative operator along the normal to the hypersurfaces as 

do = —{d t -£$). 

The vacuum evolution equations have the form 



(3) 
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d Q Kij = -^d k d k ij + ^d k d( ij)k + ^g ab d {l d labU) - Q + a^j g ab d {l d j)ab 

+ (g ab C a{ll)b +i 9lJ C - e-QdidjifP) + Tlij , (4) 
dodkij = -2d k K i:j +rig k{i Cj) + X9ij C k + T^kij , (5) 
where the constraint variables 

C = h g ab d k (d a hk — dkab) + (Hamiltonian constraint) 
Cj = d a K a j — g ah djK a b + Ti-j, (momentum constraint) 
Ckij = d k ij - d k g tj , (definition of d kij ) 

Ci k ij = d[id k }i 3 ■ ("closedness" of d kij ) 

have been added to the right-hand side of the evolution equations with some free parametcres CiIjViX- Here, 
d k = g kl di and bj, d k are defined as the traces of d k 



i j 



bj = g" l d ki] , d k =g v d ki j. 

The "remainder" terms 1Z, IZi, IZij and lZ k ij are homogeneous polynomials of degree 2 in dkij, diQ and (this 
will have direct consequences on the linear stability, as discussed below), where the coefficients may depend on gij. 
Finally, the Lie derivative of the symbols dkij is 

£ fjd k ij 

The evolution system has characteristic speeds {0, ±1, ±-\/Ai, ±\/A2, ±\/A3}, where 

Ai = 2a, 

A 2 = 1 + X- 2^ + C)»? + 7(2 -V + *X), 
A 3 = ^X+^(l-C)r?-i(l+2a)(77 + 3 X ). 

It should be noted that these are the characteristic speeds with respect to the do operator. This means that A — 0, 1 
correspond to propagation along the normal to the hypersurfaces or the light cone, respectively. The characteristic 
speeds with respect to the dt operator are obtained from these after the transformation 

/in Nfi + iVX , 

where n l is the direction of the corresponding characterisitic mode. 

The conditions under which the system is completely ill posed (CIP), weakly hyperbolic (WH) or strongly hyperbolic 
(SH) (see next subsection for these definitions) were found in KST. The system is CIP if any of the above speeds is 
complex, while it is WH if Xj > for j — 1,2,3 but one of the conditions (|^, |7]) below is violated. Finally, the system 
is SH provided that 

Aj > Q, fori = 1,2, 3, (6) 
A 3 = ^(3A! + 1) ifAi = A 2 . (7) 

For example, if the parameters (C, 7, rj, \) are ai l zero the dynamics is equivalent to the ADM equations written 
in first order form with fixed densitized lapse and fixed shift (which are WH). If a = as well, then the system is 
equivalent to the ADM equations with fixed lapse and shift (which are also WH). 

In KST seven extra parameters are introduced and used to make changes of variables in and dkij ■ When 
performing this change of variables the constraint Ckij = is also used in order to trade spatial derivatives of the 
three-metric for d k ij- Thus, the equations with the new variables have different solutions off the constraint surface. 
However, one can see that at the linear level this change of variables does not involve the addition of constraints. 

B. Linear stability 

Here we study the linear stability of the KST system considering perturbations of Minkowski spacetime written in 
Cartesian coordinates. That is, we assume that the background quantities are 

= S tj , Kij = 0, d klj =0, N { — 0, Q = 0. 



4 



In fact, since both Q and N l can be freely specified, we will first fix them to their background values for simplicity (the 
more general case is analyzed at the end of this section.) Because the non-principal part of the evolution equations 
for Kij and dkij depends quadratically on Ky, dkij and the derivatives of Q, only the principal part remains when 
we linearize these evolution equations. More precisely, the linearized equations have the following structure: 

9ij ^Kij , 

Kii = (Ad)ij, (8) 
dkij = (BK)kij , 

where A and B are spatial, first order, differential linear operators with constant coefficients. Furthermore, pertur- 
bations of the three-metric do not appear in the equations for and dkij- As a consequence, it is sufficient to 
consider only these equations; after having solved them, the three-metric is obtained through a time integration. So, 
the relevant linear evolution equations have the form 

3 

u = Y, Aj 9jU, (9) 
i=i 

where u is a "vector" formed by the components of and dkij (u has, thus, 24 independent components); and A 3 
are 24 x 24 matrices. 

Since the matrices A 3 have constant coefficients and we want to make an analysis in the absence of boundaries, 
we consider a three-torus with periodic boundary conditions as a domain and analyze equation (j^) through Fourier 
expansion. That is, we write 

u (t, x) = J2 u(t, %y t3 . 

k 

The solution to equation (^|) is, then 

u(t,x) = ^e 1P ^V^u(0,fc) , 

k 

where P(k), called the symbol, is defined by 

3 

P(k) :=J2 Ajk i- 
j=i 

There are three different possibilities for the behavior of e lP ^ k ' t : 

1. The system is SH: P has real eigenvalues and is diagonalizable, the solution simply oscillates in time. Never- 
theless, we should keep in mind that the three- metric can grow linearly in time if k = 0, see Eqs. (||). It is not 
difficult to see that this mode corresponds to an infinitesimal coordinate transformation. 

2. The system is only WH: P has real eigenvalues but is not diagonalizable. Besides the zero-frequency growth in 
the metric, the Jordan blocks of dimension n + 1 in P allow for growth in u that goes as (kt) n , where k = \k\. 

3. The system is CIP: P has at least one complex eigenvalue. Besides the previous growing modes, there can be 
exponential — frequency dependent — growth, i.e. u w exp (constant x kt). 

To summarize, the linear stability properties of Minkowski spacetime depend only on the hyperbolicity of the 
system, since the non-principal terms are zero (except for the evolution equation for , but this equation decouples 
to linear order.) Stability advantages of the choice of parameters empirically found by KST in their analysis of a 
single black hole spacetime cannot be displayed in the example we have considered. 

Finally, we study the effect of having a non-vanishing linearized (maybe densitized) lapse and shift: In this case, 
there are two extra terms in the evolution equation for Kij and dkij '■ 



Kj = (Ad)ij - didjQ, 

d ki] = (BK)kij + 2d k d( i N j ). 



(10) 
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Since we assume that Q and N l are prescribed, the solutions to these equations will have the form of the sum of a 
homogeneous solution and a particular solution. It is not difficult to find a particular solution, because we expect that 
a variation of Q and N % induces an (infinitesimal) coordinate transformation i 11 h i^ 1 + with respect to which, 
in this case, 

9ij i-> 9ij + 29(iXj) , 

Kij h-> Ku-didjX*, (II) 
dkij >— ► dkij + 2dkd(iXj) . 

In fact, if we make the ansatz Kij = —didjf, dkij — 2dkd(i£j), we get a particular solution to the system ( flO| ) if 

/ = 2ad l Z l + Q, 
4 = dif + N;. 

After a Fourier decomposition, this system has the form 

v(t, k) = iM(k)v(t, k) + q(t, k), 

where v and q are formed by the Fourier amplitudes of /, £j and Q, Nj, respectively (J — 1 . . . 3), and 



M(jfe) 



( 


2crfci 


fcl 





k 2 





\h 












By choosing the homogeneous solution to satisfy the appropriate initial data, it is enough to consider a particular 
solution with v(t = 0, k) = 0, 

t 

v(t, k) = e iM W [ e~ lM ^ s q( Sl k)ds. 



If a > (as must be the case if the system is SH), the matrix M is diagonalizable, with real eigenvalues 0, ±y 2a\k\ 2 . 
Since the norm of exp(iM(k)t) can be bounded by a constant that does not depend on t nor on k (for example, if 
2<T = 1, the matrix M(k) is symmetric and exp (iM(k)t) is unitary), this implies that v can grow at most linearly in 
time if q is uniformly bounded in time. 

If v grows at most linearly in time, the same will hold for the main variables <7y, and dkij, see (|TT]) . Again, we 
do not see any dependence on the parameters found by KST in this particular example. 



IV. NUMERICAL EXPERIMENTS 

In H we show analytically that the iterated Crank-Nicholson (ICN) method with any number of iterations and the 
second-order Runge-Kutta methods do not yield convergent discretization schemes for systems of equations that are 
not strongly hyperbolic, even if dissipation is added. Here we exemplify this through the numerical evolution of the 
KST system in a simple situation. We evolve (£|^|,^]) in one dimension. That is, all quantities in those equations are 
assumed to depend only on t and x. 

The sense in which stability is used in this section is that of the Lax equivalence theorem, i.e. a scheme is said to 
be stable if the numerical solution u£ at time t = nAt satisfies 

IKII</(*)II«°II 

for all initial data u® and small enough At and Ax. Here the index n corresponds to the time step and k to the spatial 
mesh point. It is in this sense that stability plus consistency is equivalent, through Lax's theorem, to convergence 
(convergence with respect to the same discrete norm in which the scheme is stable and consistent). Note that in this 
context consistency means that the discrete equation approaches the continuum one in the limit of zero grid spacing. 
In what follows we use the discrete norm 

ikii 2 = E«) 2a ^ 
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We will use second order Runge-Kutta (RK) for the time evolution. The spatial derivatives were discretized with 
centered differences, adding explicit lower order numerical dissipation as well, the amount of it being arbitrary. Most 
of the simulations presented here will be for RK, though the results are very similar for ICN (see the end of the 
section). In order to isolate our simulations from effects coming from boundary conditions, we wil l ch oose, as in 

and of Ref. 



Ill 



Section III, periodic boundary conditions. Also, even though the analytic stability analysis of Section 
|| is linear, here we will present nonlinear simulations, but of situations that are close to flat spacetime. In all cases 
we will use the same initial data, same numerical scheme, and change only the continuum formulation of Einstein's 
equations. 

Writing the evolution equations as 

u = Au' + l-o. , 

where "l.o." stands for lower order terms and time integration done with RK corresponds to 

1 + C ( 1 + ■ 



ul +1 = 



'k - 



(12) 



where 



C = ^S - Id-XS 4 + At x l.o. , 
with A = At/ Ax the Courant factor, a the dissipation parameter and / the identity matrix, and 

S U k = Uk+l - Mfc-l , 

5 4 u k = Uk+2 - 4tifc+i + 6u fc - Au k -i + u k - 2 ■ 



(13) 
(14) 



Our strategy for comparing numerical stability with the level of hyperbolicity is the following: we fix all parameters 
of the KST system except one, and change that single parameter in such a way that one gets a SH, WH or CIP 
system. 

The metric that we evolve, giving the corresponding initial data, is 



ds 2 = e fl+f2 (-dt 2 + dx 2 ) + dy 2 + dz 2 , 



(15) 



where f\ = f\(t + x) and f% = fofy—x) are arbitrary functions. This is an exact solution of Einstein's equations that is 
pure gauge, i.e. it is flat. In order to illustrate our point it is enough to make a simple choice for these functions, such 
as fi = A sin (x + t) and fa — with x € [— tt, it]; furthermore, for the simulations shown below we choose A — 0.01. 



A. Densitizing the lapse 



A necessary condition in order to get well posedness within the KST system is to choose a positive densitization 
for the lapse, a > 0. Furthermore, one has to take a — 1/2 if one wants physical characteristic speeds. It is worth 
noticing that with a = 1/2 and the metric (|l5|), the perturbation in the densitization of the lapse Q is zero; thus, as 



shown in Section 



III 



even off the constraint surface the solutions to SH formulations should only oscillate in time, 
except for possible zero-frequency linear growth. 

In a similar way, choosing a = results in a WH system, and a < in a CIP one. 

Densitizing the lapse is not enough for well posedness. One also has to add the constraints to the evolution 
equations. For definiteness, in this subsection we will show simulations using ( = —1,7 = 0,r) = 4, \ = (which 
corresponds to the EC formulation without making the change of variables), but the main conclusions do not depend 
on this particular choice. 



The SH case (a = 1/2) 



In order to understand some of the features of the SH numerical results that we will present we start analyzing, 
from a discrete point of view, the same numerical scheme when applied to the model equation 



u = u' 



(16) 
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Inserting the discrete Fourier mode 

into the difference scheme (|l2|), one obtains the discrete dispersion relation 



• A/ — 1 + ( i\sin{f3Ax) - 165 sin 4 ^ 



1 + sin(/3Aa;) - 85 sin 4 ^ 



p(/3Ax) 



(17) 



(18) 



where S = a\. Equation ( |18| ) should be seen as a relation between uj and (3. The discrete symbol, p((3Ax) tells us 
how the different Fourier modes are damped and dispersed (see, e.g. M, chapter 7). Writing to = a + ib, we have that 



e~ bAt = \p(/3Ax)\, 



(amplification factor) 



a 
1 



(3At 



arg(p(/3Ax)) 



(numerical speed) 



(19) 
(20) 



Choosing, for example, A = 1/2 and S = 0.01, one can see [|2| that the scheme ([L2|) for Eq.(|l6[) is stable and, therefore, 
convergent. From now on we will use, unless otherwise stated, these values. Figure [l] shows the magnitude of the 
amplification factor (|n]) and the numerical speed ( po|) for this choice of parameters. 

The continuum equation (^|) is neither dissipative nor dispersive. Ideally, one would like the difference scheme to 
have the same properties. For small values of A the initial data essentially consists of only one Fourier mode, Asinx, 
which corresponds to (3 — ±1. If one uses, say, 60 gridpoints and a spatial domain that extends from — tt to 7r, then 
one has that j3Ax — tt/30, for which Eq.(|l9|) predicts a damping of about 0.0032% per crossing time. Similarly, from 

Amplification factor Numerical speed 

X = 1/2, S = 0.01 (SH model equation) I = 1/2, S = 0.01 (SH model equation) 




j3A.x 




FIG. 1: Plot of the amplification factor and the numerical speed associated with difference scheme (|l^). If the product (3 Ax is 
small enough, then the damping and the error in the propagation speed of the correspondent Fourier mode are very small. 



Eq.(p0|) one can see that for this resolution and this initial data one should have a wave loss (i.e. the numerical and 
exact solution have a phase difference of 2ir) at approximately t — 4, 600, which corresponds to 750 crossing times. 
Figure |^ show both the numerical and exact solution for the g xx component of the metric at different times. The 
predicted speed difference is there and is in very good agreement with the analytical calculation based on the model 
equation. Also, the profile of the numerical solution maintains its shape for a very long time. Figure ^| shows the L2 
norm of the errors for this component at the same resolution. The errors come mostly from the speed difference, and 
the maxima and minima correspond to a phase difference between the analytical and numerical solution of odd and 
even multiples of n, respectively. Indeed, if one computes the error in the metric due to a phase difference of n, one 
gets 

^sin (t+*)/10Q _ e sin (t+^/lOO^ 2 ^ ^ = Q (335 

which is in remarkable agreement as well with the maxima in figure |^. This speed difference is a rather usual 
numerical feature in solutions of finite difference schemes approximating hyperbolic equations. Usually this difference 
is smaller for bigger Courant factors S , as can be seen from figure for the particular discretization here used. One 
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FIG. 2: Exact and numerical solutions for the SH case, at different times, showing the speed difference. This run was done with 
60 grid points, A = f/2, a — 0.02 and a domain of length 2n. With this resolution, the code 'loses' one wave after, roughly, 750 
crossing times. 



could try to push this factor to its stability limit, or use a different scheme that minimizes this effect. But here the 
point that we want to emphasize is that the errors in the SH simulation are second order convergent ones, as shown 
in figure |[ That figure shows the Li norm of the error in one component of the metric, g xx , for resolutions ranging 
from 120 to 480 gridpoints, up to 1, 000 crossing times (in order to evolve up to 1, 000 crossing times without loosing 
a wave, one has to use more than 60 gridpoints). One can see that the errors grow linearly in time. One could be 



tempted to think that this is due to the zero- frequency linear mode predicted in Section III. This is not the case 



in this simulation this mode is not excited (though it could be, in a more general evolution) and the error is caused 
by the numerical speed difference, as discussed above. Indeed, by performing a Fourier decomposition in space of 
the numerical solution one explicitly sees that the amplitudes are roughly constant in time (for non zero frequency 
components, this is exactly what the linearized analysis at the continuum predicted). Figure [| shows some of these 
amplitudes for this simulation (to be contrasted later with the WH and CIP cases). 



2. The WH case (a = 0) 

Figure [7] shows plots of the errors associated with evolutions performed with the same initial data, dissipation 
and Courant factor as above, except that now we densitize the lapse according to a — ("exact lapse"). For two 
fixed resolutions the coarsest one gives smaller errors after a while, and the time at which this occurs decreases as 
one increases resolution. This indicates that the difference scheme is not convergent. Note that one could be easily 
misleaded to think that the code is convergent if one did not evolve the system for long enough time, or without 
high enough resolution. For example, if one performs two runs, with 120 and 240 gridpoints, one has to evolve until, 
roughly, 150 crossing times, in order to notice the lack of convergence. To put these numbers in context, suppose one 
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Errors in the metric and extrinsic curvature 
Densitizing the lapse (SH case, 60 gridpoints) 

0.050! , , , , , , , , , , , , , , , , r 




time 



FIG. 3: 1/2 norm of the errors for the metric and extrinsic curvature in the SH case. These errors are mostly from numerical 
speed difference, as described in the text. 

Numerical speed 
S = 0.01 , Ax = tc/30, p = 1 (SH model equation) 




I 

FIG. 4: Dependence of the numerical speed on the Courant factor. Increasing the Courant factor, decreases the error in the 
propagation speed. 



had a similar situation in a 3D black hole evolution. To give some usual numbers, suppose the singularity is excised, 
with the inner boundary at, say, r = M , and the outer boundary is at 20M (which is quite a modest value if one 
wants to extract waveforms). In this case 120 and 240 gridpoints correspond to grid spacings of, approximately, M/5 
and M/10, respectively (usual values as well in some simulations). If one had to evolve up to 150 crossing times in 
order to notice the lack of convergence, that would correspond to t = 3, 000M, which is several times more than what 
present 3D evolutions last. Of course, the situation presented in this simple example need not appear in exactly the 
same way in an evolution of a different spacetime, or with a different discretization; in fact, in the next subsection we 
show an example where the instability becomes obvious sooner. Also, there are some ways of noticing in advance that 
the code is not converging. Namely, it seems that the numerical solution has the expected power law growth that the 
continuum linearized analysis predicts until all of a sudden an exponential growth appears. But if one looks at the 
Fourier components of the numerical solution, one finds that there are non-zero components growing exponentially 
from the very beginning, starting at the order of truncation error (see figure ||) . 

One might expect that, since for WH systems the frequency-dependent growth at the continuum is a power law 
one, it is possible to get convergence by adjusting the dissipation. In || we show that even though certain amount 
of dissipation might help, the code is never convergent and, indeed, adding too much dissipation violates the von 
Neumann condition, which leads to a much more severe numerical instability. We have systematically done numerical 
experiments changing the value of a without being able to stabilize the simulations (more details are given below), 
verifying, thus, the discrete predictions. 
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FIG. 5: L/2 norm of the errors for the metric. 

Fourier components of the metric 
Densitizing the lapse (SH case, 60 gridpoints) 
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FIG. 6: Fourier components of the numerical metric for to = 0, 4, 8. 



3. The CIP case (a = -1/2) 

Figure ^ shows the error in the metric, for different resolutions. As in the WH case, the errors originate mostly 
from the non-zero frequencies (i.e. the ones that typically grow in an unstable numerical scheme). But now they grow 
more than 10 orders of magnitude in much less than one crossing time and it is quite obvious that the code is not 
converging. This is so because in the CIP case the instability grows exponentially with the number of gridpoints. This 
can be seen performing a discrete analysis for the single ill posed equation in ID, v t = iv x . One gets that the symbol 
p(PAx) is real and cannot be bounded by 1 in magnitude, making the difference scheme unstable (independently of 
resolution). If one changed to characteristic variables exactly this equation would appear in ID as a subset of the 
system that we are considering, so this model equation is, in the linear approximation, part of the system evolved in 
this subsection. The amplification factor for this model equation is plotted in figure [H]. 

B. Adding the constraints 

Now we fix all parameters except the one that corresponds to the addition of the Hamiltonian constraint, and 
choose three values that will render SH, WH and CIP systems. The results are very similar to those of the previous 
section, so we will present them in a less detailed way, except that now we will discuss the role of explicit dissipation. 
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FIG. 7: Z/2 norms of the errors for the metric. 



Fourier components of the metric 
Densitizing the lapse (WH case, 60 gridpoints) 
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FIG. 8: Fourier components of the numerical metric for u = 0, 4, 8. Some of the components grow exponentially from the very 
beginning. 



1. The SH case 



We start by considering the one-parametric KST subfamily of strongly hyperbolic formulations of Eq.(2.39) of KST: 

<r=\, C = -^(23 + 20 X ), 7 = ~ V=l- (21) 

This family has, for any value of x, characteristic speeds or ±1 (see Section pll]) , for the simulations here shown wc 
chose x = 1 (and A = 1/2, a = 0.02, as in the SH simulations previously discussed). Figure [l^ shows the analog of 
figure ^. As before, we do have convergence. The Fourier components of the solution remain constant in time, and 
the error comes mostly from the numerical speed difference. 



2. The WH case 



Here we also use the parameters given by ( pl| ) with \ = 1, except that now we choose 7 = —32/21. With this choice, 
Ai = A3 = 1, but A2 = 0; therefore, as summarized in Section III, the system is WH. As in the previous subsection, the 
difference scheme may appear to be convergent, but in fact it is not, see Figure |l^. Frequency-dependent exponential 
growth starting at very small values accumulates and causes the instability. We have exhaustively experimented with 
different values for the numerical dissipation without being able to stabilize the code. Next we show some plots to 
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FIG. 9: Li norm of the errors for the metric. 



illustrate this. In figure |l4| we plot one of the Fourier components of the numerical metric. We increase the dissipation 
parameter a, starting from a = 0.02, and double it each time, while keeping the resolution and Courant factor 
fixed. At the beginning the rate of the exponential growth becomes smaller when one increases cr, though it is never 
completely suppressed, which causes the code to be non convergent. However, when one reaches the value a = 0.32 
the instability is even worse than adding less dissipation, and the same thing happens if one keeps on increasing a 
beyond 0.32. So we next narrow the interval in which the dissipation is fine tuned, we start at a = 0.24, and increase 
at intervals of 0.01. We find the same result, at a > 0.25 there is already too much dissipation and the situation is 
worse. Fine tuning even more, we change a in intervals of 0.001 around 0.250, but it is also found that for a > 0.250 
the effect of more dissipation is adverse, as also shown in figure [l4]. 

The fact that beyond a = 0.250 the situation becomes worse is in perfect agreement with the discrete analysis of 
. There we show that a necessary condition for the von Neumann condition to be satisfied is a A < 1/8. Here the 
upper limit of 1/8 corresponds to, precisely, a = 1/4. Exceeding this value results in a violation of the von Neumann 
condition; as explained in when this happens there is a numerical instability that grows explonentially with the 
number of gridpoints (i.e. as in the CIP case), much faster than when the von Neumann condition is satisfied (in 
which case the growth goes as a power of the gridpoints). 

Finally, it is worthwhile to point out that we have also tried with smaller Courant factors, using, in particular, 
values often used in numerical relativity, like A = 0.20 and A = 0.25, without ever being able to get a completely 
convergent simulation. 



The CIP case 



Finally here also use the parameters ( pl| ) with % = 1, but now we take 7 = —79/42, which implies A2 = —1 and, 
thus, the system is CIP. The results are as expected. There is exponential, frequency-dependent growth that makes 
the numerical scheme unstable, see Figure UB. 



C. Other simulations 



Performing simulations with ICN instead of RK yields similar results, as predicted in ||. Figure |l6| shows, for 
example, evolutions changing the densitization of the lapse, as in the first subsection, but using ICN with two 
iterations (counting this number as in |]). This is the minimum number of iterations that yields a stable scheme 
for well posed equations but, as shown in it is unstable for WH systems. The same values of Courant factor and 
dissipation as above were used in these runs. We have also tried with other values of Courant factor and dissipation 
parameter, finding similar results. We were able to confirm the lack of convergence predicted in Q in every WH 
or CIP formulation we used, including the ADM equations rewritten as first order equations in time space. Lack of 
convergence with a 3D code, using the ADM equations written as second order in space and ICN, for the same initial 
data used here, has also been confirmed [fTo||. 
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FIG. 10: The rate of growth of the Fourier components increases as one increases the resolution. 



V. DISCUSSION 



We have shown that a linearized analysis of the KST equations implies that flat space-time written in Cartesian 
coordinates is a stable solution of the equations if the parameters are chosen in such a way that the system is strongly 
hyperbolic. No further restrictions on the parameters are placed by this analysis. We have also integrated numerically 
the KST equations and shown that the system cannot be made stable (in the sense of Lax) and therefore convergent if 
the parameters are not chosen in such a way that the system is strongly hyperbolic. No amount of artificial viscosity 
was able to fix the problem. 

The conclusions of Section III are, to some extent, similar to those of Alcubierre et al, who find that, at the linearized 
level, the advantages of the Baumgarte-Shapiro-Shibata-Nakamura |ll[] with respect to the ADM formulation come 
from the fact that the first one is 'more' hyperbolic. However, since it is not 'completely' hyperbolic, an ill posed (zero 
speed, in the notation of |l2| ) mode is still present in BSSN. The conformal traceless (CT) decomposition is therefore 
introduced as a way of decoupling that mode. However, one could, instead, get rid of that mode by using a different 
choice of lapse. More explicitly, in the analysis of |l2| exact lapse is used, but it is not difficult to see that choosing, 
for example, densitized lapse with, say, a = 1/2, gets rid of the ill posed mode (indeed, certain version of BSSN with 
densitized lapse has been recently used in a 3D evolution of a single black hole |lj|). Moreover, as shown in Jl4) , 
the result is much stronger: even at the nonlinear level appropriately writing BSSN with densitized lapse results in a 
reduction (from first to second order in space) of certain strongly or symmetric hyperbolic formulations. 

The lack of convergence in evolutions of WH or CIP systems with schemes that are convergent for well posed 
equations seems to have been overlooked in the past. If the von Neumann condition is satisfied, one could be easily 
misled to think that the scheme is convergent, especially if coarse resolutions are used, the evolution is short, or few 
frequencies are present in the initial data. Here we have shown that this lack of convergence does appear in concrete 
simulations, even in very simple ones. These numerical experiments, together with the theorem of reference should 
be enough evidence to cast serious doubts on any simulation performed with evolution equations that are weakly 
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FIG. 11: Amplification factor associated with the difference scheme (Jl2|) approximating the ill posed equation v t = iv : 
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FIG. 12: L2 norm of the errors for the metric. 



hyperbolic that are not accompanied by very detailed convergence studies. 

Summarizing, we do have analytical tools at our disposal to constrain schemes and predict to a certain extent their 
behavior under numerical evolution. In fact, the examples shown suggest that in the case of linearized equations, the 
analytical tools are complete: one can a priori tell if a code will work or not. In the non-linear case further tools will 
need to be developed to achieve the same status with respect to predicting the performance of a code before running 
it. 
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FIG. 13: Z/2 norm of the errors for the metric. The simulation is stopped once the determinant of the spatial metric becomes 
zero. 
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FIG. 14: We were not able to make the code stable by fine tuning the dissipation parameter. 



Errors in the metric 
Adding the Hamiltonian constraint (CIP case) 
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FIG. 15: 1/2 norm of the errors for the metric. 




FIG. 16: L2 norm of the errors for the metric using ICN. 



